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Abstract 

The detection of B-mode polarization in the CMB is one of the most important outstanding tests 
of inflationary cosmology. One of the necessary steps for extracting polarization information in the 
CMB is reducing contamination from so-called "ambiguous modes" on a masked sky which contain 
leakage from the larger E-mode signal. This can be achieved by utilising derivative operators on the 
real-space Stokes Q and U parameters. This paper presents an algorithm and a software package to 
perform this procedure on the nearly full sky i.e., with projects such as the Planck Surveyor and future 
satellites in mind; in particular, the package can perform finite differences on masked, irregular grids 
and is applied to a semi-regular spherical pixellization, the HEALPix grid. The formalism reduces to 
the known finite-difference solutions in the case of a regular grid. We quantify full-sky improvements 
on the possible bounds on the CMB B-mode signal. We find that in the specific case of E and B-mode 
separation, there exists a "pole problem" in our formalism which produces signal contamination at very 
low multipoles I. Several solutions to the "pole problem" are presented; one proposed solution facilitates 
a calculation of a general Gaussian quadrature scheme, which finds application in calculating accurate 
harmonic coefficients on the HEALPix sphere. Nevertheless, on a masked sphere the software represents 
a considerable reduction in B-mode noise from limited sky coverage. 



1. INTRODUCTION 



The detection of primordial gravitational waves in the Cosmic Microwave Background 
(CMB) would be the "smoking gun" of inflationary cosmology[l, 2]; the existence of such 
gravitational waves in the inflationary scenario implies the observation of primordial tensor 
pertubations in the microwave sky in the form of curl-type ("S-mode") polarization. When 
using harmonic methods on a masked sky, signal-mixing with the divergence-type ("£'-mode") 
polarization occurs[3], as the E- and B-mode power contributions from the masked regions 
are unknown (and hence the contributions from masked regions and their boundaries are also 
known as "ambiguous modes") and the signal decomposition into polarization modes is in- 
herently non-local. Since the cosmological 5-modes are realistically expected to be at least an 
order of magnitude smaller than the £^-modes[4], such signal-mixing could potentially either 
drown out any true B-mode signal or masquerade as a false positive detection. Instead, one 
may perform such full-sky masked extractions by utilising real-space derivative operators on 
local analogues of the E- and S-modes via finite dijferences[5] and thus avoid sampling in the 
masked regions at all. The finite-differencing method is a standard method for performing 
numerical derivatives, where derivatives are estimated by attaching numerical weights to the 
input data at a set of pixels surrounding each individual pixel. 

The most popular pixellization for CMB analysis is the semi-regular spherical HEALPix grid [6]; 
standard HEALPix methods for extracting E- and S-modes are based on harmonic techniques 
and thus suffer from signal-mixing when a mask is present[7]. Performing numerical deriva- 
tives on irregular grids, over a range of geometries, is a problem which is both general and 
non-trivial. While other, more regular, grids such as GLESP[8] exist for CMB analysis, they do 
not hold the same computational advantages as HEALPix. 

We present the MasQU (Masked Stokes Q, U analysis) software which can perform such deriva- 
tives via finite-differencing. The weights are not derived for particular grid types and derivative 
orders as in the standard approaches, since a completely general method to calculate them al- 
gorithmically is developed. The algorithm is then tested on a Cartesian grid, and on the sphere 
(as supplied by the HEALPix scheme). The sphere is known for being problematic due to the 
poles; it is non-trivial to define the position differences for a pole-crossing pixel sample since 
the coordinate pair {9, 0) is not unique at either pole. Even after this is dealt with, in the par- 
ticular case of extracting E and S-modes there exists an accuracy problem related to the pixel 
positions in the HEALPix polar cap. 
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We find that, compared with standard harmonic techniques our method reduces leakage by 
at least an order of magnitude, the reduction increasing with multipole scale. Compared with 
the flat-sky approximation often used, the full-sky formalism is a vast improvement on large 
angular scales. 

The paper proceeds as follows: in section II the CMB polarization formalism is summarized 
and we develop the finite-difference weight-generation scheme, quantifying pixellization leak- 
age. Section III then presents an implementation on the HEALPix spherical grid — including 
discussions of some of the pitfalls in attempting to perform accurate derivatives at the pole. The 
performance of the software on the HEALPix grid in the presence of both noise and masking 
is analysed in Section IV, in comparison with the standard harmonic methods. We also make 
mock calculations for EBEX-type experiments. Finally, we conclude with some comments. 

II. ANALYZING CMB POLARIZATION ON THE FULL SKY 

Standard inflationary scenarios in cosmology predict the existence of tensor fluctuations in 
the modem Universe. Tensor fluctuations from inflation are "frozen-in" to the CMB on su- 
perhorizon scales; this means that if primordial S-modes exist they will be most prominent at 
small values for the multipole /, corresponding to large angular scales. In order to fully probe 
large angular scales we need a full-sky analysis, something the Planck Surveyor[9] will provide 
data for. In fact, whilst £^-modes have been detected since 2002 (DASI[10]), the cosmological 
S-modes have remained elusive — the current state of detection from combined Wilkinson 
Microwave Anisotropy Probe (WMAP) analysis giving a tensor-to-scalar ratio r < 0.24 at 95% 
confidence[ll]. In the following section we describe our formalism for decomposing the mea- 
sured CMB polarization into parameters useful for cosmology. 

A. CMB Polarization Formalism 

The Stokes parameters Q and U, which describe the linear polarization of light, are de- 
pendent on the choice of reference basis but can be reparameterized in terms of two basis- 
independent polarization modes known as the "E" and "B" modes[12] by analogy with electro- 
magnetism, since any tensor field of rank greater than zero can be decomposed into "gradient"- 
type and "curl"-type parts. The relation between the two sets of parameters on the spherical 
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surface is given by 

Q{n) = J] (oL hYimin) + -2Yim{n)] + iat hYUn) - -2Yim{n)]) 

^ (1) 
U(n) ^--J2{aL hYimin) + -2YUn)] + iaf^ [2YUn) - -2Yim(n)]) , 

Im 

where the gYim are spin-weighted spherical harmonics[13] with integer spin weight s — reduc- 
ing to the standard spherical harmonics Yi^ for s = — n is the chosen coordinate basis and 
the are the s = harmonic coefficients of X = {E, B}. Since one can decompose vector and 
tensor fields into curl and divergence parts, and cosmological vector fields decay exponentially 
in the inflationary scenario, the decomposition of tensor perturbations corresponds uniquely to 
the E- and S-modes. In the real universe of course, the situation is confused by the presence of 
lensing, reionization and other astrophysical phenomena restricting the detection of primordial 
tensor mode phenomena to low-/ multipoles. 

To facilitate a real-space calculation, we can construct scalar and pseudo-scalar fields corre- 
sponding to the E and B-modes respectively 

^(^) - E ^Jl^.^'irnYiUn) h{n) = ^ Ji^.^irnYUn) (2) 

and then relate the bi-Laplacians of these fields to the Stokes parameters (where we have 
dropped the n for convenience, while 5" = 5 • • • 5 and 8" = 8 • • • 8): 

V^e = -^[B^Q + iU) + Q\Q - lU)] V% = '-[B\Q + iU) - Q\Q - iU)] (3) 

where the 8 and § terms are respectively the Newman-Penrose spin-weight-raising and spin- 
weight-lowering operators 

— —{de -\- i CSC 9dtt> — s cot 9) 8s = — {d^ — i esc 9d^ -\- s cot 9) (4) 

and we have appended an s to the standard notation which allows easier spin-weight-counting 
of the entity operated on; the operation of 8s on a quantity with spin-weight s returns a quantity 
with spin- weight s and correspondingly the operation of 8 returns a quantity with spin- 
weight s — 1. Rearranging equation (3) gives 

V^e = -D^^Q - D-,U V% = D^.Q - D+,U, (5) 
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where the derivative operators are 

5_iS_2 + SiS2 

— —A — CSC ya.^^ + 6 cot voq + oee 

(6) 



-0^2 = — ^- — ^ ^— = -2 - csc^ 6*9,^0 + 3 cot Qde + 



-D=p2 = ^ ^^"^ ^"^^^^ = 2 esc 6'(cot Bd^ + a^,^) 

which can be used to extract e and h from the underlying Stokes polarization maps, with the 
referring to the combination of spin-weights of the (S, 8) operators which define the operation. 
The "inverse" relations for forming Stokes fields from the underlying e, h fields are 

Q = -Die + D-h U = -D+h - D'e (7) 

where 

Do = ^ ^0 = ^. ■ (8) 

In the flat-sky small-angle approximation, since the spin-weight s of a function gf is defined by 
its transformation on rotation about the pole via 

./ ^ e-^/, (9) 

all spin-dependence of the operators vanishes 

5 = -(a, + idy) 8 = -(a, - idy) (10) 

yielding the approximation 

D^2 ~^ Dq ~ dxx — dyy ~^ Dq "^dxy (11) 

For these operators, the relations 

[D+,D^]=0 (flat sky) 
D+^Do - D-2D+ = (full sky) 

are satisfied; the non-commutativity of their discrete analogues is a measure of signal leak- 
age from pixellization. Since it can be shown that the harmonic-space contribution of the bi- 
Laplacian is 

^ (13) 

then relating the power spectra of the scalar and pseudo-scalar fields to that of the E- and 

5-modes is trivially 

_ (/ + 2)! g y4j, _ (/ + 2)! ^g 
^' - (I^^' ^' - (I^^' ■ ^^^^ 



(12) 



5 



B. Computing Derivatives by Finite-Differencing 



A derivative of a function f{x) with respect to coordinate x with order m at a given pixel 
labelled i can be computed as the sum of weighted values of the function at a surrounding 

sample of pixels j[5]: 

pixels 

j 

where w is the weight matrix. The finite-difference method involves estimating derivatives by 
taking weighted differences of function values at selected pixels (see Appendix A for further 
details); a general method[14] for computing the weights on an arbitrary Id grid requires, for 
the derivatives at a pixel i, solving the linear system 



/ AO 



A" 



A? 



\ 



\^ A-^ • • • A-^ J \ J \{n- l)lSm,n-i ) 



( 



(16) 



with a Kronecker delta for separating the m*^ derivative such that only a single entry on 
the right-hand side is nonzero for the given derivative, and A^^ = Xj — a position difference 
term. Various tests of this method can be found in Fig. 1. 

The geometric array in equation (16) is the differenced "Vandermonde" matrix[15]. The Van- 
dermonde matrix is related to the Lagrange basis interpolating polynomial 



Li{x) 



n 



(17) 



which can be used to interpolate a function / sampled at n points by 



(18) 



In the Vandermonde formalism this polynomial can be identified with the determinants of ar- 
rays determined by the pixel sample taken about i: 



Det[t;,] 
~ Det[v] 



(19) 



where the subscripts indicate that the Vandermonde matrix defined at pixel j has the column i 
replaced by a column of undetermined positions x, i.e., for the 1-dimensional case 




jpt fed diff scheme vatidii] 




3pt central diff scheme varied i] 




-(.! 0,0 O.i 
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4pt 1 -hii diff scheme, vaned ,V| 




-1,0 -0,5 0,0 0,5 1,0 



-1,0 -0,i OO 0,5 1,0 



ipt bwd diff scheme, varied ,ri 
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FIG. 1. Accuracy 5x = Xgxact— a^caic of 2-, 3-, and 4-point finite-difference schemes on regular and irregular 
Id grids (positions on the grid are xq < xi < X2 < x^, assuming existence), for a first-order derivative. 
We use the test function /(x) = x^, since there are fewer datapoints available in these schemes than 
can specify the function exactly. The differences blow up only as the coordinate variation produces 
degeneracy in the datapoints. 



X 



V2 



(20) 



X 



n—l ^n~l 



X. 



n-1 
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in order to find the interpolated function value at a given pixel i (not necessarily part of the pre- 
sampled pixel distribution), with x corresponding to the positions Xj in equation (17), placed 
in column i in Vi. The Vandermonde array is then the 1-dimensional case of a more general 
geometric array, related to interpolation over d dimensions. The specification of the geometric 
array is dictated by the dimensionality of the space being interpolated over and by inspection 
of the geometry (and thus possible derivatives) of the pixel sample one is working with; for 
example, a square grid of 4 x 4 points in {x, y) allows a polynomial specified by powers of 
up to x^, and combinations thereof; the geometric array would then include combinations 
of powers of x, y up to such a limit. For the purpose of performing derivatives, we replace 

Xi y ^ij ^i' 



C. Estimators 



For cosmology, it is necessary to calculate the errors on these terms in harmonic space such 
that one can define the accuracy of the power spectra estimator. For example on a flat Id space 
we can compute the estimator for the transform of a 1^* -order derivative, 

k k 

where fk is the Fourier power in F for mode k. On a discrete grid there is instead the 1-point 
forward difference equation 

k 

where the ~ over the derivative operator denotes a discrete derivative, hence the finite- 
differenced spectrum Vk always underestimates the true spectrum Vk'- 

Vk = sinc2(A;A)Pfe. (23) 

In the following, we define an n-point stencil as a collection of n pixels over which the finite- 
difference calculation in equation (15) for a given point on a grid is made. The full flat-sky 
operators on a regular grid with 3-point stencils are then given by 

2 

A2 



' 2 (2^) 

k 



For more general Id derivatives, we have to start from 

~d,F = /fc^^'" E ^j^''""''' (25) 

k j 

and induce simplifications at this point using grid symmetries. There is then no simple sinu- 
soidal statement for general irregular grids but the estimator is not difficult to calculate. The 
1^* -derivative spectral estimator is then proportional to 

Vk oc [^Wje^^^i^i'^ (26) 

where the full irregular flat sky discrete operators are found using the general solution given 
by equation (A28) in appendix A 



^ij = - 



PIJ 



(27) 



(for focal pixel i, number of points n, order m) with the operator 9a defined by equation (A29), 
which can then be used to compute the leakage. 

For the full-sky discrete estimator, an algebraic solution to the 2d geometric array is highly 
involved. Instead, one can compute (by brute force) the numerical power contributions for each 
derivative term via the Wigner 3/m-symbol[16] by treating the Wj elements as part of {n + 1)^ 
separate fields (for the case of 2-dimensional square grids): 

{dn)lm = ^Jl^f\^i)yi',rn'm,,^Y,Un,) (28) 

j,l',Tn' "^^^ i 

where the Qj^i refers to the rotations between the focal pixel i's position on the sphere and that 
of the given neighbouring pixel j. It is then possible to calculate the accuracy of the power 
contribution by measuring the bi-Laplacian (the "signal" generated by the derivatives) and the 
commutator (the "residual" generated by the derivatives): 

Sig ^ D^,D^ + D-,D^ = Res ^ D+,D^ - D'.D^ (29) 

since equation (5) can be written in pixel space as 

A . , ,\ . / (30) 

V% = (p^^D^ + 1^+2^0 ) b + Err(6) + [D^^D, - D^,D+j e, 

where Err represents remaining error in the calculation from the summation of the "signal" and 
"residual" parts. Whilst we can take a simple approximation by assuming leakage into e from 



b is negligible the full calculation is clearly coupled. A calculation for the residual estimator, 
using the scheme described in the next section, can be seen in Fig. 2. This calculation is then 
equivalent to defining the pixelization leakage in terms of a convolution operator 

{dn)lm — y^ll'mm'Yl'm' (31) 
I'm' 

where 

yVll'mm'= J (^W,{n)^Yi,m'imimmd^, (32) 

with which one can deconvolve the pixellization leakage. In fact, it is will be useful to derive 
the more general convolution operator 

,..w„w = /5*l..r,w(n)l,i',„(n)<in (33) 



if one wishes to perform the exact calculation in the Stokes tensor frame (see Appendix B). 
Returning from these more general considerations, the residual in the Stokes tensor frame can 
be approximated by calculating the operators involved in equation (3) 

S_i5_2 = dee — csc^ QdA,A, + 2i esc QdaA, + 3 cos B esc Qde + 2i cos Q csc^ QdA, — 2 

(34) 

5i52 = dee — csc^ ^5^^ — 2i esc Qde^ + 3 cos Q esc Qde — 2i cos 9 esc^ — 2 

which are spin-weighted in correspondence with operation on the Stokes tensor instead of the 
individual Q and U fields, and their discrete analogues, given by 

5_i5_2/(0) = 5i52/(^^) + R-^fm ^id2f{ft) = 8182/(0) + R-fm, (35) 

where the residual operators i?^ are the quantities we will approximate. By taking the lowest 
order errors in the Taylor series derivation of the first and second order derivatives in Id (see 
Appendix A) 

^2 ^ 12 

'—^^SJ' + ^9,,f' (36) 

Ji+1 ' Ji-1 Ji-1 Ji+1 Cl £.i , ^ 



the operators are 



R^^AK -de* ± i—des4, + —des ) - -^^t>' ^ -^de^^ T 2i cos csc^ Od^^ \ (37) 



10 



and by substitution into 



Im 



Im 



W% [^""-2^;^ - R-2Ylm\ + ^ J] C [R'^-2Ylm + R-2Ylm\ 



(38) 



Im 



Im 



the residual approximation is calculated. For the case of HEALPix, one has A^(6',0) 
Ao{9, (f)) I sin Q, yielding the contribution at each point in Q\ 



E 

_ Im 



1+4)! 
(/-4)! 

^(/ + 4)! 



«L I 1 + 



6sin^ 



1^ 



E 

Im 



''Im 



V 6 sm e* / 



(/ + 4)! 
(Z-4)! 

'(/ + 4)! 



1 + 



6sin^ 
1 

6sin^ 



1^ 



Im 



(39) 



Alternatively to utilizing equation (28), one might simply reconstruct the Q and U fields solely 
from the sum of the underlying e-field harmonic coefficients and perform the real-space com- 
mutator derivatives. We will also see in the next section that there is a consistency criterion for 
calculating S-modes on the real-space HEALPix sphere which can help inform us as to whether 
we have calculated a signal which is dominantly numerical noise. 
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FIG. 2. Left: Fractional signal (black lines) and residual leakage (red lines) of the finite-difference scheme 
on the HEALPix sphere (as described in sec.III) at O2; (Thick line, dashed line) corresponds to A'side = 
(8, 16). There is a characteristic high point at the very \ow-l scale, followed by an immediate dip in 
leakage. The first phenomenon is related to the pole problem in our formalism on the HEALPix sphere; 
the signal starts to peak again at approximately the Nyquist frequency for the map (~ 2A'side). Right: 
Residuals from individual point Ci maps generated from = dw where I' specifies the point Z-range. 
In this case, we have /' = 10 for maps with A'side = 8, 16 (thick, dashed lines). The leakage has a 
steep downward gradient across / with the leakage contribution to each multipole oscillating rapidly, 
the oscillation rate scaling with map resolution. 



III. THE HEALPIX SPHERE 

HEALPix[6], a Hierarchical Equal Area isoLatitude Pixelisation scheme for the sphere, is the 
most widely used software for construction and analysis of full-sky CMB maps. The lowest- 
resolution partitioning is into 12 equal area pixels. Each pixel is assigned a unique identification 
number (in either a "nested" or "ringed" numbering scheme), and is surrounded by 8 other 
pixels, except in the polar cap where some of the pixels are surrounded by 6 or 7. The map 
resolution is specified by the parameter 

iVside = 2^-<^- iVorder G |Z| (40) 

with each map of 12N^^^^ pixels composed of 4A^'side — 1 isolatitude rings; at each level higher in 
resolution, the pixels are subdivided into 4 equal area pixels in the higher-resolution map. The 
rings immediately by either pole consist of 4 pixels (independent of resolution), increasing by 
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4 pixels per ring for each ring increment toward the equator, up to a maximum 4A^side pixels in 
the equatorial rings. 

The HEALPix sphere is an interesting case study, since it is a semi-regular distribution on a 
coordinate system with a pathology: the multi-valuedness of at the poles — hence differences 
m.{6,(j)) across the pole are ill-defined. It does not matter if we do not directly sample the pole 
point itself — merely crossing the pole with a differencing stencil is enough to complicate the 
calculations, since the interpolating polynomial covers the whole region. Similarly, we must be 
careful at the — 0/27r boundary. However, at such a boundary here we can simply reassign 
the differences to the smallest of the two possible paths across the sphere between the points 
we are differencing. 

To test the software, the derivatives of harmonic functions on the sphere are computed, since 
primordial cosmological point signals are not expected. This implementation uses LAPACK at 
double precision and the truncated SVD technique[17, 18] for ill-conditioned matrices (in the 
limit of large stencils and irregular geometries). The geometry for each pixel stencil taken is 
given by a square-geometry differenced geometric array 



/ 1 



1 \ 

A0„ 
A^„A0„ 

A^r'A^r' / 



(41) 



where we have chosen to work in the {6, 0) basis for computational convenience (one does not 
expect the Stokes fields to necessarily be a polynomial in this basis). This choice is determined 
by issues arising from analysis at the pole; whilst the analytic operators we are calculating are 
necessarily covariant, discretization complicates matters at the pole. Consider a function f{9, (p) 
which is smooth in {6, 0); at the pole there is a coordinate singularity in {6, 0), forbidding the 
use of finite differences in this region. A change to Cartesian coordinates {x, y) would swallow 
the coordinate pathology; however it can be shown that g{x, y), the Cartesian-space description 
of f{6,(j)), will exhibit a functional pathology at the origin {i.e., g {0,0) may be a delta function or 
similar). This phenomenon is general for coordinate transforms which remove the polar coor- 
dinate singularity. Since we are approximating our function with a polynomial, any functional 
discontinuity or pathology will be poorly modelled by such a finite-difference scheme. The 
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utilised alternative is choosing the pixel stencils in {9, 0) such that they do not cross the pole at 
all, but instead progress toward an outer-difference scheme around the pole. 

In the particular case of extracting E- and fi-modes, the pole also presents other numerical 
problems; one way to deal with these problems is by "rotated oversampling" around the pole. 
This will be discussed briefly after the next subsection. 

A. Structure of the Algorithm 

The algorithm is as follows: 

• Construct an approximately square stencil of nearest-neighbour pixels, one for each pixel. 
This is achieved by taking the unique set of nearest neighbour pixels for the focal pixel, 
and then repeating recursively for the neighbour pixels until the specified stencil radius 
is satisfied. If the array cannot be filled with corresponding pixel numbers (such as in 
the case that a pixel has less than 8 neighbouring pixel), then the remaining elements are 
assigned an identification value of -1. Note that since the HEALPix nearest-neighbours- 
finding routine always calls the neighbour pixels in the same geographic order, then pixels 
with identical surrounding geometries will necessarily have identical stencil arrays. In 
other words, the code is structured in such a way that the symmetries of the HEALPix 
grid are preserved also in the stack of pixel stencils taken. Our notation for the stencils 
is as follows: a stencil of order 0„ contains at most (existence and masking-dependent) 
precisely (n + 1)^ pixels. Since we wish to bias toward central differences, each n is even, 
corresponding to a radius of n/2 pixels about the central focus pixel. The stencil order n 
is initially the same for all pixels. 

• Perform "re-mapping" for pixels pi surrounded by one or more masked pixels — search 
amongst the stencils of the surrounding pixels pj for the "optimum" stencil (the stencil 
with the most available pixels and closest to the central focus). A weight system of 



is applied. Pixels without enough information in their stencil to perform the requisite 
derivatives (three available pixels with a unique position 9, and the same for ) are dis- 
carded. 




If pixel is missing or in a cut region 

1 Otherwise 
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• Analyse the mapped difference geometries to single out a smaller number of unique sten- 
cil geometries — these characterize a unique set of stencils which provide all the necessary 
weights for calculation across the sphere, saving computation time. In the unmasked case, 
there are also symmetries between the north and south polar regions and also the quarters 
of each hemisphere which can be taken advantage of to cut down computation time. 

• Find solutions to the linear equation 

v'w = 5, (42) 

whose Id analogue is equation (16), corresponding to each of the derivatives required — 
each set of weights is unique for a given stencil geometry and derivative order (except in 
trivial cases). The geometric array is by default that for a regular square array. In the case 
that there is not a full square stencil available, we remove rows (starting from the lowest 
row /largest difference powers) and the corresponding columns of the geometric array 
such that the array remains square. It is probably more rigorous to remove rows using 
geometric considerations instead (i.e., an analysis of which derivatives can be calculated 
from the stencil), but this has yet to be implemented and does not significantly effect the 
rest of the analysis since only the derivatives up to second order in {9, 0) are calculated. 
Furthermore, the additional time cost for such a procedure might not be a good trade-off. 
Solutions to the linear equations are found using LAPACK and the QR decomposition or 
SVD technique, depending on the pixel sample. 

• Repeat until all the calculations for each pixel are finished, and form the bi-Laplacians. 
Compute the power spectra of the resulting e and b maps by removing the power con- 
tributed by the bi-Laplacian operator. 

• In the presence of a mask, one should apodize the signal since masking redistributes 
signal power. An apodization subroutine is available with the software, whilst there is an 
optimal method for CMB studies due to Smith & Zaldarriaga[19]; we leave the analysis 
of this to a later paper. 

• Output the derivative maps, the power spectra, and the weights used (these can be recy- 
cled once calculated for the first time) in FITS format. 

The above methods are very easily tweaked in the software for a given geometric scheme (be 
it for alternative pixellization schemes on such as GLESP, or non-spherical schemes) and a 
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given desired derivative. Higher-dimensional generalizations have not been implemented but 
would be trivial. Testing the software for accuracy follows the process of taking pre-specified 
E- and B-mode harmonic coefficients, calculating the full-sky Q and U maps from these coeffi- 
cients, operating on the sky maps to produce the bi-Laplacians, and then converting the power 
spectra of the e and b fields to that of the E- and 5-modes for comparison with the original 
input spectra: 

dQ,dU 

where we also compare the full scalar field maps V^e, with the sum of the original harmonic 
coefficients as in equation (2). The accuracy of the software, as operated on a test harmonic 
function, is shown in Fig. 3. 

The SVD technique separates a matrix A into 

A = (43) 

(where S is the diagonal array of singular values) and is known for yielding optimal linear 
solutions in the presence of near-singular matrices. In the near-singular case, the smaller singu- 
lar values can be dominated by round-off error leading to dramatic errors in the solution. For 
large geometric arrays, the bottom-row elements are most likely to suffer from round-off error 
since they are large powers of small numbers. This means the singular value array must be 
truncated below some numerical threshold in order to yield reasonable numerical results, and 
replace the corresponding inverted elements in the inverted singular value array with zero; 
truncation error will contribute to any inaccuracy of the calculation. Since different resolutions 
correspond to a scaling of geometries on the HEALPix sphere, the threshold should be defined 
by the ratio Ejj/En. Optimal truncation may depend on some non-machine aspects: array size 
and geometry. Whilst this problem has not been solved generally, empirically-derived trunca- 
tion thresholds for the unmasked sphere as-a-whole (i.e., individual stencil geometries have not 
been studied) which maximize the accuracy of the calculations have been implemented for the 
first few stencil orders. Masking has not been fully SVD-optimised (we use the same thresh- 
olds as for unmasked maps) since this study requires generalization to any number of pixels in 
a range of distributions. 

Finally, we consider time complexity. While it can be seen that a first-time implementation of 
the code is out-performed by the standard HEALPix methods in terms of speed, one needs 
calculate the weights only once. For a survey like Planck, with up to 5 x 10^ pixels per map, 
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FIG. 3. O2 and O4 accuracy (across rows) for derivatives d0,deg,d(f),d^if,,de(f, (down columns) of 

T 

/ 1.5 X IQ-^ 0.1 4.5 X 10"^ 2.6 x lO^^ 3.1 x 10"^ \ 
^22(6*, </>)/ with absolute maximum values of ~ . 

\ 1.5 X 10^2 Q 8 15 X ^ l X 9 Q ^ iQ-2 J 

Since the error in the polar cap is much larger than the equatorial region error for all 
maps bar 9^, the values in the equatorial regions in these maps has been limited to 

T 

( ±10-6 ±2 X 10"*^ ±2 X 10-6 ±2 X 10-6 ±10-6 \ 

. This shows that the accuracy improves with 

y ±10-3 ±10-3 n/a ±10-3 ±2 x 10-^ / 

stencil size, and is worst at the pole (due to a combination outer-differencing error (for de and dee), and 
increased differencing values for the A terms). However, it is not reasonable to simply increase the sten- 
cil order with the aim of achieving some threshold accuracy, for three reasons: 1) The time complexity 
for linear solutions goes as 0{it'), where n is the number of pixels in the matrix; 2) We are limited to 
a maximum stencil size of (2A''side ± 1)^- (3) The q;5ample (/, m) given has a relatively small magnitude 
pole problem error for harmonic functions (See Fig. 5). Thus the convergence rate with stencil size is too 



the maximum HEALPix resolution is Ai'sidc = 2048. For the corresponding O2 calculation, we 
used an interpolating function to extrapolate a timescale of 1440 seconds on a low-end machine 
(Fig. 4) since for a fast performance the procedure requires arrays of size ~ h{n + l)^A^pix, which 
exceeds the machine's RAM limit. 




2 4 6 8 2 4 6 8 

Log2[Nside]-2 Log2[Nside]-2 



FIG. 4. Left: Time complexity for a maskless (black lines) calculation across Nsi^^- The key is {O2, O4, 
06) = (Thick line, long-dashed line, short-dashed line). Due to RAM limitations we use interpolation 
and extrapolation for the Planck resolution timescale, to2,Pianck = 1440s, to4,Pianck = 9024s, toe.Pianck = 
39946s. Masked calculations are only very marginally slower (and the time costs display some small 
mask-dependency). Right: By comparison, HEALPix calculations with default parameters. The grey line 
is the zeroth order calculation, with the black line for 4 iterations of the HEALPix map2alm subroutine. 
For a Planck resolution map the timescales are 534.7s and 1823.36s respectively. 

B. The "Pole Problem" and Rotated Sampling 

For polar pixels, the algorithm is amended as follows: 

• Define polar pixel "re-mappings"; for each pole-crossing pixel (referred to as Ppoi), we 
search among its stencil pixels for the set of pixels whose stencils do not cross the pole. We 
then select from this set the subset of pixels that satisfies the requirements that minimize 
their distance on the sphere from ppoi in order that the new calculation scheme is as close 
to a central-difference scheme as possible, and from this subset choose the pixel which is 
closest to the pole (referred to as Pnpoi)- We then reassign the stencil of Ppoi from a central- 
difference in the Ppoi stencil basis to an outer-difference in the chosen pnpo\ basis. 
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The full array of pole problems for schemes not crossing the pole are as follows. Firstly, we 
can see from our 1 -dimensional tests (Fig. 1) that while the forward- and backward-difference 
schemes are of the same magnitude error as the central schemes, they generally perform slightly 
worse. The same is true for comparisons of higher- to lower-order derivatives. This is because 
higher-order derivatives require more sampling pixels, and hence a larger basic stencil size, 
than lower orders. We can also see this in the HEALPix scheme for derivatives of the s = 
spherical harmonics (Fig. 3), hence the higher errors at the poles. We refer to this as "outer- 
differencing" error. 

Secondly, the A0 values between sampled pixels increase toward the pole due to the lower pixel 
sampling per HEALPix ring toward the pole. At its most extreme, the pixel rings unmediately 
surrounding either pole have only 4 pixels each, with a separation 7r/2. Since the errors in our 
differencing scheme depend on powers of A0, then we expect a lower accuracy around the 
pole. This we refer to as the "polar A0" problem. It is apparent in all cases that (i) Second-order 
derivatives perform slightly worse than Ist-order derivatives; (ii) dee performs worst. This is 
no surprise given that the polar cap pixels are only central in 0; (iii) large-magnitude functional 
variation in about the pole seems to correlate with large polar errors (Fig. 5). 

We can look first at ameliorating the outer-difference and polar A0 error. Fig. 6 includes 
calculations when the stencil size is increased, but only around the pole. 
A third, more drastic problem comes from the construction of the terms we wish to calculate. 
In the continuum limit, these are 



V^e = -(-2 - csc^ edss + 3 cot Ode + dee)Q - 2 esc ^(cot Ods + des)U 

(44) 

V% = 2 CSC ^(cot ed^ + de4>)Q - (-2 - csc^ Od^^ + 3 cot Ode + dee)U. 



This is quite a delicate combination due to a number of esc 9 and cot 9 terms, which clearly 
blow up as we approach the pole. The contributor of the largest error is then the csc^ 98^^ term, 
which will blow up any errors in the discrete approximation to d^^. This is referred to as the 
"blow-up" problem, and accentuates the first two problems. 
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FIG. 5. Top 2 rows: O4 accuracy maps of 6-fields for a range of unit-a^„ sources; the maps are symmet- 

/ 5.9 X 10-'^ 0.6 133 137 \ 
ric in errors, with maximum absolute values of . Bottom 2 rows: Q maps 

\ 0.8 0.3 3.8 611 J 
corresponding to the sources given. Whilst the absolute maximum values are of order ~ 1 for each Q 

map, the polar errors peak for high m = 2 and m = 3, and the magnitude of the polar errors at fixed 

tn increases with /. (/, m) = (32, 32) maps have been included to show how error can propagate at the 

equator, due to variation in cj) greater than that modellable by the stencil. 

The combination of these effects leads to dramatic problems at the pole for realistic (CMBFAST- 
generated[20]) sky maps: in Fig. 7 we can see the operation on realistic i?-mode-free maps. 
Even for a large stencil the error swamps the e signal; for a non-zero b signal the problem is 
then more drastic due to the b signal's small magnitude. 

Taking derivatives using larger stencils is not always advantageous. Instead, a method for 
dealing with the pole is by rotated sampling; we rotate the sphere by some small amount 
S(p < A0 and perform derivatives on a "doubled" stencil. Fig. 6 shows that the error is greatly 
reduced in this manner. There is though the caveat that it is necessary to sample the harmonic 
coefficients accurately in order to effectively remove the pole problem. The polar problem 
ensures that taking "discrete" doubled sampling (calculating over a stencil which includes the 
next-neighbour-in-0's stencil) is not nearly effective enough. Of course, for the rotated sam- 
pling case we can always perform the rotation more than once in order to improve the accuracy. 
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In the case of discrete doubling, we are limited to using A^pix-in-ring/2 adjacent stencils, so as 
not to have a pole-crossing calculation. 
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FIG. 6. O2 6-fields: Various oversampling techniques; only the polar values are of interest here. From 
top to bottom: (i) Original results for the source af^ = 6i2Sm2, af^ = 0, (ii) with an O4 stencil about the 
pole, (iii) Scj) = 0.01 HEALPix-reconstructed rotated doubled sampling, (iv) 64> = 0.01 analytic rotated 
doubled sampling, (v) "discrete" doubled sampling. The absolute maximum values for each map are 
1.20 X 10^ 116 3.03 X 10"^ 104 186 ) • The "discrete" doubled sampling mode performs better than 
the standard O2 but falls short of the rotated sampling; with an accurate map reconstruction method 
there is nothing in principle to forbid an n-tupled rotation sampling, which would bring the pole error 
down quickly. 
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FIG. 7. Top row: O4 V^b map generated from CMBFAST Q and U maps, with tensors. Bottom row: O4 
map generated from CMBFAST Q and U maps, no tensors. The map resolution is A'^^ide = 128, and all 
the other parameters are the defaults from the LAMBDA[21] online tool page. Columns, left-to-right: No 
ring removal, 1st north & south polar rings removed, equatorial region only. In the no-tensors scenario, 
at the power spectrum level the pole problem could provide a false-positive detection; at the level of 
the map, we can distinguish these by eye — inconsistency between polar and equatorial regions in a 
given map, as seen in the middle and right-hand images, is a result of the differencing error dominating 
the calculated map. This then amounts to a consistency criterion by which to check for a false positive 
detection. The irregular geometries outside the equator induce a larger error than at the equatorial 
region. For the equator, the cut-off length is defined by the ring at which the stencil geometries become 
irregular, and is hence O^-dependent. 



C. The "Pole Problem" and Accurate Reconstruction of the Harmonic Coefficients 

We know^ that the construction of our bi-Laplacians is highly sensitive at the pole to error 
in the underlying derivatives, and that rotated sampling can solve this problem. In order to 
produce an adequate rotated map, one must accurately calculate the harmonic coefficients aim 
of the map, apply the Wigner rotation functions, and then sum over the new coefficients. 
The HEALPix method of reconstructing the harmonic coefficients of a scalar field is an iterative 
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procedure [22], manifested in the map2alm subroutine: one starts with a zeroth-order estimator 



(45) 



and resums the coefficients to form a map T^^\ The next step is to take the difference map 
^2^(0) _ 2" _ 2^(0) Qj^^ compute the zeroth order harmonic coefficients of 5T^^^ in the same manner 
as the zeroth order of T, then iterate and sum over coefficients to form the n^'^-order approxi- 
mation: 



-1 Np 



(46) 



3=0 i 



The optimal aim sampling scale for map reconstruction is /max = 2A^sidc/ with an optimal number 
of iterations being 3 according to the HEALPix software recommendations. While this is a quick 
and reasonable approximation, a numerical analysis finds that the convergence limit for the 
iterations is not suitable for the rotated sampling method. 
We wish to perform the integral 



aim = J F{n)Yii{n)dn 



(47) 



For a HEALPix grid, the points in 6 are described by[6] 

4 2k 



cos 9 — 
cos 9 



3 

1 - 



3-/Vside 



North equatorial belt, A^side < ^ < 2A^side 
North polar cap, 1 < k < A^side 



(48) 



(and correspondingly in the southern hemisphere). Although this is globally irregular, the 
HEALPix 9 points may be split into regular (equatorial) and irregular (polar cap) parts. 
For the regular grid part one might choose a Newton-Cotes [17] method. Since the equatorial 
region is equally spaced, we can use the composite trapezoid rule for the poles and sum it with 
a single ^equator-point Newton-Cotes (NC) approximation at the equator. 
If an integral is approximated by 



J fdx = E ^ifidx, 



(49) 



where the nodes are equally spaced, then it can be shown that the NC weights can be found by 
solving 



xl 



X' 
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(50) 
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Table I shows the results for this method on a regular grid, and Table II a comparison with 
Gaussian methods. 



TABLE I. Ntunerical test of the NC scheme; from Mathematica we have j\f{x)dx ^ 1.71125; A-^ = h—a. 



No. of nodes 


X 


NC scheme 


result 


2 


-1,1 


^(/o + /i) 


1.21306 


3 


-1,0,1 


^(/0 + 4/l + /2) 


2.60653 


4 


-1,-1/3,1/3,1 


^(/o + 3/1 + 3/2 + 4/3) 


2.1771 


5 


-1,-1/2,1/2,1 


^(7/0 + 32/1 + 12/2 + 32/3 + 7/4) 


1.71047 


6 


-1,-3/5,-1/5,1/5,3/5,1 


^(19/o + 75/1 + 50/2 + 50/3 + 75/4 + 19/5) 


1.71082 


7 


-1,-4/6,-2/6,0,2/6,4/6,1 


^(41/o + 2I6/1 + 27/2 + 272/3 + 27/4 + 2I6/5 + 41/6) 


1.71128 



Gaussian quadrature on the HEALPix sphere 



Since we know how the regular and irregular HEALPix 6 grid parts are constructed, one can 
then attempt a mixed NC and Gaussian scheme for the regular and irregular parts respectively, 
or even a Gaussian scheme for the whole sphere. Gaussian quadrature (in particular Gauss- 
Legendre quadrature) is used in GLESP to yield fast and highly accurate harmonic coefficients 
without recourse to iterative methods. The HEALPix points do not follow the abscissas of the 
Gauss-Legendre scheme, so we need to derive the general Gaussian scheme first. Starting from 
the integral[23] 

F{x)W{x)dx f=i J2 ^jfi^j) (51) 
and using an interpolating polynomial Lj as given in equation (17), one finds 



Wi 



-dx 



such that the relevant integral to solve is 



dx. 



A general solution to the above integral with W{x) = 1 over m points is 



m—l 



Im = ln(a; - Xj) Yl{xj -Xi) + ^ 



x^ 



i=l 



z=l 



a=0 



a 



combinations 



•^yo ' ' ' -^ya ) 



(52) 



(53) 



(54) 
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giving for example a 3-point solution 

/3 = y + Y - ^ Xi^ + X ^1X2 + X1X3 + X2X3 + a;| - Xj ^ Xi^ (55) 

where the log term disappears since we can only realistically sample at the pixel positions i. 
The performance of this method is better than the NC method, and almost competitive with 
Gauss-Legendre (see Table II) — although it can be shown that our general Gaussian scheme 
reduces to Gauss-Legendre at those nodes. 



TABLE II. 4-pt results siuximary for various methods 



Method 


value 


Exact 


1.71125 


Gauss-Legendre 


1.71122 


Newton-Cotes (regular) 


2.1771 


General Gauss Method (regular) 


1.7222 


General Gauss Method (irregular) 


1.67456 



There does however remain an important problem with this scheme: the number of com- 
binations to find for, say, an A^side = 32 HEALPix map for a mixed Gauss and NC (Gaus- 
sian at the polar caps) and full Gaussian scheme respectively are Yfa=i 31! [(31 — q;)!q;!]~^ and 
^126^ 127![(127 — o;)!a!]~^ which have prohibitively large time complexities. 
A final alternative is to rotate the underlying coordinate system hy 59 = 7r/2 which rotates the 
pole to the equator of a new coordinate system. This allows us to calculate the rotated scalar 
and pseudo-scalar field pole pixels using a central-difference scheme whilst avoiding the pole 
issues of the esc ^-type terms. Rotating the resulting calculations back to the pole results in 
ringing from the pixel boundary, which can be tempered by calculating a larger region of ro- 
tated pixels since ringing dies down further from a discontinuous boundary. Performing such 
a method essentially doubles the computation time of the full field calculations, but is recom- 
mended if the polar pixels are needed. 

In the absence of a both quick and accurate method for performing rotated sampling, one is re- 
duced to simply removing the offending pixel rings; for most calculations this is recommended 
since the polar pixel region is negligible for large (i.e., Planck-type) maps. A comparison with 
various oversampling techniques is seen in Fig. 6. 
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IV. LEAKAGE 



Our analysis so far has been limited to an ideal unmasked sphere; however the real sky is 
obscured by the galactic plane amongst other foreground sources and will have various contri- 
butions of noise to it. We use this section to note the performance of the software on a maskless 
sky against the standard methods (extending the results of section II.C), i.e., pure pixellization 
and finite-difference error, and then show the efficacy of the software in the presence of mask- 
ing and noise. Of course, whilst testing against standard methods shows that the proposed 
method works as required, the real test of the software will be against alternative proposals [19, 
24, 25, 26, 27, 28, 29, 30] for clean subtraction of the i?-mode leakage. 

Figure 8 compares MasQU calculations on the full sky to similar HEALPix calculations; on 
the full sky the harmonic space separation is superior to calculations using the first few stencil 
sizes; it being the case that with large sky masking the MasQU method is superior (see next sub- 
section), there should exist an intermediate masked area, for a given MasQU calculation order 
and map resolution, wherein the difference between the HEALPix leakage and MasQU error 
is minimized. We leave this interesting calculation to a future publication, since for realistic 
experiments the masking volume will likely exceed such an equilibrium value. 
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FIG. 8. Maskless calculations, with HEALPix on the left and O2 MasQU on the right. For the HEALPix 
calculations, the black line is for iterations while the grey is for 3 iterations of map2alm. For the MasQU 
calculations, the grey line is for a standard calculation, with the pole removed for the black line. The 
corresponding e signal has an rms ~ 2 x 10^ ^K; the red lines are full-sky B-mode spectra from tensor 
modes corresponding to r = 10^^, 10^^, 10^^ (left-hand diagram) and r = 10^^, 10^'^, 10^^ (right-hand 
diagram). 
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Pole removal (Fig. 9) has the effect of lowering power across scales, with the decrement 
most noticeable at low /. For completeness, we also show the difference between operating in 
the full-sky and flat-sky formalisms; it should not be surprising that the errors in the flat-sky 
approximation are at their largest at low /. 
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FIG. 9. Left: O2 and O4 (black, grey) maskless S-modes from the r = map, with the pole removed. 
Right: O2 calculations with the flat-sky and full-sky operators (black, grey). In both plots, the red lines 
are B-mode spectra from tensor modes corresponding to r = 10^^, 10^'^, 10"^. 



A. Masking 

In the case of masking, analysis was performed using 3 basic types of mask (Fig. 10): equa- 
torial, polar and random masks. 

In the masked case differences are found as would be expected; the pole problem particularly 
affects the power at around / = 2 by boosting it significantly. Focussing on the no-tensors maps, 
at multipoles higher than / ~ 10 MasQU performs considerably better than the raw pseudo-C/ 
calculations for masked HEALPix schemes, by about 1 to 3 orders of magnitude; in most mod- 
els gravitational lensing rather than primordial modes dominate the foreground polarization 
from / ~ 150, meaning we have a large /-range where MasQU is advantageous for calculating 
5-modes. The smoothness of the masked calculations is in contrast to that of underlying func- 
tional discontinuities. The effect of a shelf discontinuity itself is shown in Fig. 12. As expected 
(since the underlying approximation to the signal is an interpolating polynomial), the software 
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FIG. 10. Rows, top-to-bottom: Masking schemes, V^e maps and maps, at Oq and A'sidc = 32 for 
^Fm ~ (^hn ~ ^- source function was chosen for the smallness of the pole problem errors, but 

nonetheless the masking errors are smaller or of magnitude that contributed by the pole. 

actually performs worse with a larger stencil when in the presence of a discontinuity. For the 
CMB, discontinuities will mostly be contibuted by point sources on the sky; these will need to 
be masked away using a source catalogue. By extension, the software could be used as a tool to 
search for discontinuities, as discussed in [31]. 

B. Noise Performance 

A number of simple noise models (Gaussian, anisotropic uncorrelated and pixel-to-pixel cor- 
related) are also analysed. Fig. 13 presents power spectra for our realistic _B-mode-free CMB 
maps with Gaussian noise added, where variations have been made in the mean value of the 
noise (at 10%, 1% and 0.1% the mean signal value) and in the scaling of the noise (via a cut-off in 
the number of multipoles generated for the noise maps). At high magnitude and large /, since 
the noise translates to a collection of point discontinuities in real space, we expect the sum of 
the derivatives of the summed signal and noise maps to feature more point source errors. The 

29 



10 20 30 40 50 60 10 20 30 40 50 60 

1 1 




1x10"' '■■ . •'■ 

10 20 30 40 50 60 

1 

FIG. 11. S-mode Qs from our fiducial i?-mode-free maps as reconstructed by the 3-iterations HEALPix 
method (black) or an O2 MasQU calculation (grey), for equatorial, polar and random masks. In all the 
plots, the red lines are S-mode spectra from tensor modes corresponding to r = 10"^, 10^^, 10^^. 



addition of our Gaussian noise models serves to boost 5-mode power fairly consistently across 
I up to the /max value that sets the smallest scale for noise variation; a lower-l noise mode cut-off 
/max results in a drop in signal power boosting for / > /max/ but still with a significant contri- 
bution. The pixel-correlated power boost is almost indistinguishable from the pure Gaussian 
boost; in contrast, the anisotropic boost is dependent on the scaling of the direction-dependent 
noise. The MasQU method can then be seen to be very sensitive to noise; one would want a 
good understanding of the systematic and foreground noise properties in order to effectively 
purify the S-modes. 
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FIG. 12. (Logarithmic) code performance in the presence of shelf discontinuities; the i?-modes are cal- 
culated from a Q = sin 6 sin (j), U = pair of maps with a cut-off to Q = i^i the southern hemisphere 
at the equator. Without a U signal, one should have no i3-modes. Left-to-right: Original Q map, O2 
and O4 V^b maps. Notice that there is a blow-up in errors for calculations across the discontinuity, as 
expected; the error actually gets worse and effects a larger region for larger stencils, and also scales with 
the magnitude of the discontinuity. This general behaviour also follows for point discontinuities. 



C. Apodization 



Any masking scheme will ensure that the l-space power of an underlying signal is redis- 
tributed across multipoles, since there are some scales which will contribute less depending on 
the mask. Gibbs overshooting will also be a problem; boundary effects ensure that "ringing" 
occurs in any harmonic approximation. This is usually dealt with by apodizing the signal — 
applying a tapering weight to each pixel contribution to the power, such that the apodization 
function tends smoothly to zero at the boundary. 

In the E- and i?-mode case, an optimal apodization scheme in the pseudo-C/ formalism has 
been derived by Smith & Zaldarriaga[19]. It can be shown by operating 5 on a window func- 
tion representing the masking that in the analytic case passing to the scalar field is equivalent 
to the "pseudo-Crwith-counterterms" method. Specifically, 



2 J v^(/ -!)(/ + 2) v^(/- !)/(/ + !)(/ + 2) 

and equivalently for the e-term, where sW = d'^W. The optimal apodization weights for a 
minimum variance estimator are calculated by minimizing the difference between the band- 
limited pseudo-Cz estimator 

^and = Yl dmC^r^'W.d,, (57) 



31 




10 20 30 40 50 60 10 20 30 40 50 60 



1 1 

FIG. 13. Left: O2 maskless S-modes from the r = map with resolution N^^^^ = 32, with noise added to 
the underlying Q, U maps and the pole removed. The Gaussian noise models in Q and U are calculated 
from a Gaussian distribution harmonic coefficients and normalized such that it's mean is some propor- 
tion of the signal mean, i.e., {N) = a {S). Left plot: The thick black line is our noiseless model, whilst 
the scaling for the dashed black, dotted black, and dashed grey lines are a = 0.1, 0.01, 0.001 respectively. 
Right plot: The thick black line is the noiseless model, with the dashed black, dotted black, and dashed 
grey lines for Gaussian, anisotropic, and pixel-correlated noise models respectively, all with a = 0.1. In 
both plots, the red lines are B-mode spectra from tensor modes corresponding to r = 10^^, 10^'^, 10"^. 

constructed using the mask map Wi where the band-limiting is defined by a binning of multi- 
poles /, and the optimal estimator 

Oband = Yl d,iC-'C''^^''C-%d^ (58) 

where dis a data vector and C is a covariance matrix (which may be band-limited); this is then 
calculated by conjugate gradient inversion (varying the spin-weighted W components inde- 
pendently) and is only equivalent to the polarization modes in the mean of multipoles m (with 
subsequent loss of phase information). The method used in [19] is somewhere between "pure" 
(without mode-mixing) and optimal. By contrast, our real fields e and b are by construction 
pure. An application of this method to our scheme is left to a future paper. 
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D. Leakage from Realistic Surveys 



In order to get a rough idea of the leakage improvement the software can bring to real data, 
we calculate the leakage for sky coverages in a simplified model of the E and B Experiment[32] 
(EBEX) survey and the bounds that may be set on the tensor-to-scalar ratio r in these simplistic 
cases; we also perform our algorithm on a mask of similar area centred on the equator. In the 
past the signal-to-noise ratio S/N has been too low to perform differencing on the real sky, but 
such projects are designed to improve S/N, making such calculations plausible. Our fiducial 
model is much the same as the parameters set for our previous no-tensors analysis, except cal- 
culated for higher-resolution maps. 

The EBEX survey is a balloon-borne polarimeter for probing the sky with a resolution of less 
than 8 arcminutes at frequency bands centered at 150, 250, 350, and 450 GHz. The sky patch 
covered by the ~1300-detector instrument corresponds to 350 square degrees. The EBEX re- 
gion is an approximately square patch, corresponding to Figure 9 of [33]. In particular, we 
perform calculations equal roughly to the EBEX sky and ground coverages, for our fiducial B- 
mode-free model adding uncorrelated Gaussian noise in Q and U at the levels 3.2fiK and 0.9/iK 
respectively. This is performed on an iVgide = 128 resolution HEALPix map — less than the 
resolution capable by EBEX, but more than enough to capture the essential low-l polarization 
information relevant to tensor modes. We start by smoothing the signal + noise map with a 
Gaussian kernel, and calculate the E- and B-mode spectra from this (Fig. 14); the noise model 
is approximated by performing the same derivatives on a smoothed Gaussian field. On the 
HEALPix sphere, the decreasing resolution per ring will serve to increase low-/ power in the 
calculations. Since one is not using the full sky it is advantageous to rotate the survey region to 
the equator, which would ameliorate such a problem. Since the pixel distributions in the polar 
cap and equatorial regions differ significantly as in equation (48), we also perform on a similar 
area mask centred at the equator. 

Mock likelihoods for the tensor-to-scalar ratio r are also computed using CosmoMC[34]; we 
hold all cosmological parameters constant except for varying r, the scalar and tensor spectral 
indices Ug, nt and the superhorizon power of the scalar perturbations XogAg. This of course 
assumes that quantities such as the reionization optical depth are known perfectly; for a real- 
istic analysis one would have to perform the MCMC calculations over a higher-dimensional 
space that includes such parameters as variables. Thus if the noise model is well-known, then 
the MasQU method provides an excellent improvement in the mock surveys over standard 
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harmonic methods. Specifically, since real-space derivatives obey linearity, one can in principle 
approximate the noiseless map by removing the E- and 5-modes calculated by a full real-space 
noise model. Given that the noise models from projects such as WMAP are computed initially 
as sky maps, this would have the advantage of separating out point sources before smoothing 
of the map is performed. 




FIG. 14. Left column, top-to-bottom: Mock EBEX ground survey, rotated ground survey, and the ro- 
tated sky survey masks. The ground and sky surveys are imposed with a 0.9fiK and 3.2fiK smoothed 
Gaussian noise component, representing detector error. Right column, top-to-bottom: B-mode spectra 
corresponding to the survey regions in the left column, in units of (/(/ -|- l)/27r)/ui^^; the black line is a 
HEALPix calculation, the grey an O2 MasQU calculation. These results compare favourably with the 
S-mode residual in more detailed EBEX analysis [35], although the noise model we use is far more sim- 
plistic. Central column: Maximum likelihood calculations of r according to the text in section IV.D; the 
black lines are the HEALPix calculations, with the MasQU calculations in red (solid lines are the fully 
marginalized posteriors, dotted lines the relative mean likelihoods). 
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V. CONCLUSIONS 



We have presented a software package for performing derivatives on the HEALPix grid, in 
the presence an arbitrary mask. The setup given is in fact quite malleable and can be adapted 
to any pathology-free grid, for any collection of derivatives. This is particularly useful in the 
case of CMB polarization studies, where we can use derivative operators to separate out the 
E and B polarization modes and avoid mode-mixing as is problematic in standard harmonic 
methods. 

In the particular case of the operators used for CMB analysis, we found that polarization mode 
numerical operations at the pole were highly sensitive to error in the component derivatives, 
leading to very large numerical noise at the lower order differencing stencils. This was found 
to be ameliorable in any of four ways (i) throwing away pixels at the pole — while this removes 
information, there is an automatic reduction in contamination; (ii) applying larger differencing 
stencils around the pole — here we lose no information but have a slow convergence rate for 
large multipole sources; (iii) rotating the grid in theta, and then rotating back the calculated field 
values — this suffers from ringing which can be minimized by essentially doubling the entire 
spherical calculation; or (iv) "rotated" sampling — whilst we lose no information, and in theory 
have a quicker convergence at large multipoles, this requires very accurate map aim map 
reconstructions in order to be useful, something that standard HEALPix methods do not sup- 
ply 

In some senses, naive speed comparisons for the MasQU software against HEALPix harmonic 
methods are somewhat inappropriate; whilst HEALPix performs considerably faster, MasQU 
requires calculation of weights only once per masking scheme,map resolution and stencil size 
— any further calculations using the same three criteria requires only a trivial summation at 
each pixel of the recorded weights, regardless of the underlying sky map. 
Work on the software is ongoing. Further studies will involve the correct apodization of masked 
maps and a full exposition of the discontinuity-finding aspects of the formalism. A beta ver- 
sion of the MasQU software is available from the first author. Use of the MasQU package for 
producing results should be acknowledged in any forthcoming publication. Bugs or technical 
issues should be emailed to the first author. 
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Appendix A: A General Finite-Difference Scheme 

This section briefly expands on the finite-difference scheme; most of the results here are taken 
from [5, 14, 36, 37]. A derivative (of order m) of a function at a given pixel i can be computed 
as the sum of weighted values of the function at a surrounding sample of pixels j: 

pixels 

d.^f, ^ Y. (Al) 

where w is the weight matrix. For example, the canonical examples of a finite-difference scheme 
with uniform spacing can be derived starting with a 1-dimensional Taylor expansion on an 

37 



infinite regular grid with separation A 

fi±i^fi±AdJi + A 
to yield for the 1st and 2nd derivatives 



2 ^xxfi 
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Since the focal pixel in this case is the central pixel fi between fi±i, this is usually referred to as 
the "second-order central difference scheme". If the grid is finite, the first and last rows in the 
weight matrix must correspond to "forward" and "backward" difference schemes respectively, 
where the focal pixel is fi^i in the second-order case. Such schemes can be constructed using 

the same Taylor analysis as in the central difference case. 

The following more general formalism shall be constructed to approximate the derivatives at a 
single pixel; in that sense we shall utilise the weight vector wj (at fixed but unlabelled pixel i) 
instead of Wij. For the set of all pixels on a pixellated grid, the weight vector wj corresponds to 
rows of Wij. 

A general finite-difference method, for any number of regular or irregular pixel schemes can be 
derived using interpolating polynomials. In the Lagrange basis, aid polynomial interpolating 
a set of n datapoints can be written as 



f{x) ^ ^fiLi{x), 



(A5) 



where the fi are the datavalues at each point Xi and the Lagrange basis polynomial is 



Li{x) 



such that on a sample of n pixels we have 
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(A6) 
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(A7) 



where clearly Li{xj) — 5ij. It can be shown that for a given nondegenerate distribution of points, 
the Lagrange basis polynomial both exists and is unique. Interpolation problems over n pixels 
arranged on a Id grid can often be expressed using a geometric progression matrix, also called 
the "Vandermonde" matrix[15]: 



V — 



\ 



(A8) 



which is constructed by rewriting the interpolating polynomial in the monomial basis, by a 
simple rearrangement of terms in equation (A6): 



Lj = ci + CiXi H = 



(A9) 



r=l 



and then filling the rows of the matrix v as appropriate, such that 



/ ^0 
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n-l \ 
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Cl 



(AlO) 



Using this, a scheme for approximating the derivatives of the interpolating polynomial can be 
constructed; if we write equation (AlO) in the compact form vc — L, noting that the positive 
definiteness of v implies that {v~^Y — {v^)~^, we can define a unique unspecified array a which 
solves the transpose equation 

a 



corresponding to the summation 



r=l 



(All) 
(A12) 



The linear equation for these interpolation weights is then 



\ 



„n-l 



X. 



n-l 



(A13) 



By isolating a single pixel of interest, with position x, and replacing the positions of the sur- 
rounding pixels Xi with the position difference 



— X 
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(A14) 



one is led to an equation (with an unspecified array W) 



v' = WU 



(A15) 



dejfined in powers of Aj 
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0!Li 



(A16) 



which can be modified for calculating the differencing weights (where l^f* is a shorthand for 
the -nf^ derivative of Lj). We refer to the array v' as the "differenced Vandermonde" array. Let 
us now consider a more general Taylor series than equation (A2) for our polynomial function L 
in order to modify equation (A16) for calculating differencing weights: 
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(A17) 



In a matrix format, since our summation is realistically limited to the {n — Xf^ derivative this is 
none other than 

/ (±Ai)° (±Ai)"-i \ / .(0) \ ( J \ 



0! 
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V 0! 
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(±Ai] 



By inverting the matrix equation (A18), we can calculate each of the derivatives of the polyno- 
mial L. It is then clear that in order to isolate a particular derivative m, one must append the 
left-hand-side of the inverted equation with a Kronecker delta thus: 

f r \ 

\ (A19) 

y (±A„)0 • • • (±A„)"-1 ) \ Lk^n J 

where we have made use of the array (iv'^j'^ Correspondingly we modify equation (A16) 
for the same purpose: 
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in order to separate out the derivative of interest. By extension, replacing the Kronecker term 
via 5m,r Sm\\m',r wWch evaluates to 1 if either m or m' is equal to r, will calculate the weights 
required to compute the summed derivative L*^™) + L^'"''-'\ Thus the array W is related to the 
vector weights w*^™) by 



n-l 



(A21) 



3=0 

One can then solve for the weight vector by isolating the derivative polynomial of choice. For 
example, the central difference system for a regular grid, where A = 1, can be obtained from 



(A22) 
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or a backward difference equation for the second derivative via 
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(A23) 



with results in agreement with the known results. A range of standard regular finite-difference 
weights is displayed in Table III. 



TABLE III. 2-, 3-, and 4-point equidistant Ist-order difference equations. 
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The extension of this to an irregular grid merely requires the reparameterization 



(A24) 



So for the derivatives at a pixel i one solves 



\ /.„M\ / 



w. 



\ A-^ • . • A-^ J \ J \{n- l)!5^,„_i J 



(A25) 



41 



A general solution to this set-up for a derivative of order m can be determined by constructing 
the LU decomposition of the inverse of the trace of the Vandermonde matrix into lower- and 

upper-triangular arrays A and T respectively: 



' = AT 



(A26) 
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The analogous decomposition for the differenced Vandermonde matrix leads to the general 
solution for an n-point finite-difference scheme in 1 dimension 
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where i is the focal pixel and j denotes weights applied to the stencil pixels and we have defined 
the following operator: 
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" (A29) 



yielding for example 
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{Ai3 - Ail) {Ai3 - Ai2) {Ais - Aa) ( Ai4 - A^) ( A^ - A^a) ( A^ - A^g) ' 

If the determinant of an n-point Vandermonde array is labelled Detfujn, then by utilizing stan- 
dard linear algebra techniques it can be shown that 



Det[i;]„ = Yli^j - xi)Det[v]n-i. 

J=2 



(A31) 
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By iteration, the determinant of the Vandermonde array is then 

n 

Det[^;]n = {xj - Xi) 

l<i<j<n 



(A32) 



which, for a given unknown xj, with positions Xi known, is precisely the factorised interpo- 
lating polynomial. It is can then be seen that the Lagrange polynomials can be identified with 
the determinants for Vandermonde arrays determined by the pixel sample surrounding pixel j 
(where we have dropped the n): 
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where the subscripts on v indicate that the Vandermonde matrix defined at pixel j has the 
column i replaced by a column of undetermined values for x, i.e., for the 1-dimensional case. 
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since the Lagrange polynomial is merely the weighted sum of the unique polynomials at each 
point on the grid. Thus in the d-dimensional case, one would also have 



Detk 



Det[v('^)] 



(A35) 



where we have to specify the correct form of the d-dimensional geometric array. One corroUary 
of the interpolating polynomial being in the form of a geometric array is that we can imme- 
diately test for the existence of a unique polynomial (in the powers we have specified in the 
geometric array) with roots at each pixel. Such a polynomial does not exist if the geometric ar- 
ray is singular. Similarly, we can construct the new result of general solutions for the derivative 
weights in d dimensions as functions of the determinant of v': 
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where the derivative operator is 
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(A37) 



with A" = (Aei, Agj, • • • , Ae^) the vector of position differences in each of the dimensions e^, 
m" the vector enumerating the derivative orders in each dimension and P^'^^ is a multinomial 
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of A which we can construct from the (n — l)*** minors of our n x n geometric array. The 
determinant depends on the geometry of the geometric array in a very specific way — this can 
be illustrated with the use of graphs, where the 2-dimensional case is used since it is easy to 
visualize. Consider that a 2d square pixel array corresponds to the following computational 
graph: 
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where the black dots indicate interpolation zeroes, and the dashes imply continuation of the 
grid scheme. One can compute the determinant of this via the recursion 
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where there is some freedom in dimensions d > 1 in how to define the a term. Specifically, we 
can assign a graphical method to this calculation; if we append each point in the graph with a 
pair of numbers representative of the powers of (ei, 62) in the geometric array 
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then we can perform each iteration in the recursion by removing the (0, 0) point in the graph, 
and replacing each pair of numbers (a, b) with values dependent on a chosen form of the a term: 
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where the r term enumerates the r**" recursion, since each recursion removes a column from the 
geometric array in powers of (Aejrfc- The choice in equation (A39) is motivated by the desire 
that there should be a point enumerated with (0, 0) at each stage of the recursion; this point 
corresponds to a pivot column about which one can restructure the arrays without altering the 
value of the determinant. At certain stages in the recursion, the grid will feature negative values 
of (a, b), from which it is not possible to directly obtain a (0, 0) point; one then reconfigures the 
grid by calculating 

Det[v"^X = /3nDet[v'^X, (A40) 
where the array v" obeys the same recursion properties as v': 

Det[^;"(')]„ = aj)et[v"^^^ 
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The /3 factor is calculated by adding the values (a', b') to each point on the graph, corresponding 
to 
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Pn= U (AeJ7/(A,j7/. (A42) 

j=r+l 

This method then generalizes in d dimensions to form a determinant like 
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dependent on the initial array, where we have used the ~ to indicate that the solution for a 
given geometry is of the general form of (A43). Since this recursion relation is equivalent to 
calculating the determinant of a geometric array in terms of the determinants of reduced ge- 
ometric arrays, which correspond to the graph of the original geometric array but with dots 
removed, one can write this as: 





(A44) 



DetK(2)]i5 = Det 



DetK(2)]i6 = Det 

This allows one to perform a more general 2d analysis in the case where the array does not form 
a square, such as in the graph 
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The recursion terminates at the determinant vlf^ , containing elements with inverse powers of 
A. An algorithm for generating the determinant from a given graph is then: 

• Remove the Udiag diagonal terms in the d-dimensional graph defined by terms (ai, • • • , a^), 
which will contribute 



n 

i<j<ndiag 



d 



eiJik 



.1=1 



The remaining terms in the graph will have values (oi — Udiag, ■■■ , aa — f^diag)', choosing 
the shortest vertex dimension in the graph, one can then reconfigure the graph which 
will contribute 



jk 



to the determinant. 



Iterate down this vertex, removing points in the graph, to add a contribution 

n 

~n[(^eJ,fe-(AeJd 

to the determinant. 



i<j 



• Repeat until the full determinant is calculated. 

In order to complete the calculation of the general finite-difference weights, it is necessary to 
determine the polynomial P^'^^ from the minors as previously stated, using the same graphical 
method. In 2d, this term is of the form 



jj^i k<l 



or more generally of type 
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By looking at the residual R on the Taylor series truncation 
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which can be generalized to higher dimensions by utilizing the Hessian form of the Taylor 
expansion (A2) 

(A49) 



ah 



where the summation convention is assumed and dab- are the partial derivative tensors of Hes- 
sian type i.e., in 2 dimensions 
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and Nab- is a numerical factor corresponding to the factorial terms in equation (A2), it can be 
shown that the error on our general n-point d-dimensional finite-difference scheme for a square 
array is of order 

Ett ^ O (^([[A^jY/'^y (A51) 

where we recover C(A^) for 1-point radius regular central schemes in d-dimensions, or C(A^") 
for n-point radius regular central schemes. For non-pathological functions such as polynomials, 
the accuracy of the differencing scheme generally improves with the number of points used 
(Fig. 1). 

Finally, we should note that there may yet be some value in using geometric arrays featuring 
inverse powers. Whilst the Taylor series does not permit negative powers in the expansion, the 
Laurent expansion[23] of a complex function 



(A52) 



(where the a„ terms are constants and 6 is a point on the complex plane) does, suggesting that 
one can deal with a pixel at a coordinate singularity by applying the hyperbolic part of the 
interpolation at that pixel. This would, for example, correspond to the following algorithm for 
the differencing weights to approximate the Ist-order derivative on a Id grid: 
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Appendix B: General Convolution Operators 



When calculating the residual for operation of n-point discrete derivatives in different frames 
(such as the Stokes tensor frame), or even for different spin-weight fields (such as if one has a 
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vector jfield) it is useful to look at the more general convolution operator 

ss'yVll'mm' = J dn[s'Yi>m'i^)] sYlmi^)dfl. (Bl) 



We briefly expound on some of the techniques available for such a calculation, without carrying 
out an in-depth example. One can first calculate the convolution operator for the derivatives in 
the analytic case by utilising the relation between the spin-weighted harmonics and the Wigner 
D-functions 



where the "little-c?" function is 
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and the angular momentum operators in R^[16] are 

L2 = - [dl + cot Ode + csc^ 9{dl - 2 cos Od^e + d^)] 

1 (B4) 

Lz — —idcf, Lz' — —id^ L± — — ^e^*"^ (de ± i cot 9d^) 

v2 

where it should be noticed that on S^, the operator is equivalent to D^. These operators have 
the power contributions 



which allows us to construct the analytic derivatives of the spin-harmonics using the recurrence 
relations of the Wigner terms (p. 90 of [16]). As an illustrative example, the recurrence relation 
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shows that the de operator redistributes power across multipoles I. To jHnish the evaluation of 
the analytic convolution operators, it is necessary to evaluate the integral 

J s'Yvm'{^)sYiMdn. (B7) 

This can be achieved by utilising the orthogonality relations of the little-d functions 

£ dl^, {e)di^, (6) sin ede = ^5,, (B8) 

and those of the Jacobi polynomials p^^"'^^ ([38], p.806), since 
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For discrete derivatives on a regular grid with sampling points separated by a length A, the 
generalized convolution operator for derivatives in is trivial to compute. For the analysis of 
derivatives in 9, we can use a perturbation expansion in the little-rf functions, dl^^,{cos{6 + A)); 
by isolating the lowest order expansions in the sinusoidal terms in equation (B9) 
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x(l + {m' - I - m + 2r + I cos9)Acsc9) + 0{A'^) 
we see that the discrete derivatives in 9 mix power across multipoles s, I and m. 
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